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ABSTRACT: Approximate solution of the ensemble representability prob- 
lem for density operators of arbitrary order is obtained. This solution is closely 
related to the "Q condition" of A. J. Coleman. The representability conditions 
are formulated in orbital representation and are easy to use. They are tested 
numerically on the base of CI calculation of simple atomic and molecular sys- 
tems. General scheme of construction of the contraction operator right inverses 
is proposed and the explicit expression for the right inverse associated with the 
expansion operator is derived as an example. Two algorithms for direct 2-density 
matrix determination are described. 
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1 .Introduction 

In the last decade a number of Post HF methods of high precision have 
been developed ffl. It is unlikely, however, that in the nearest future these 
methods in their present form have any chance to be applied for electronic 
structure calculations of extensive molecular and crystalline systems. Instead 
much more simple DFT methods based on the work of Kohn |J are 
becoming the main tool for calculations of such systems. 

There exists an approach in density functional theory which progress 
had been restricted by serious mathematical problems connected with the 
so-called representability property of 2-electron density. After the famous 
theorem of Coleman [|J who solved the representability problem for 1-density 
operators, there appeared a number of papers with attempts to generalize 
Coleman's result to treat q-density operators (for q > 2) || ||, 0, §, |9], [1C 



These attempts, however, did not end in results of practical importance. It 
even became accepted that the problem is far too complicated to find any 
applications for electronic structure calculations. 
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An alternative attempt to use reduced density matrices for molecular 
calculations is developed in |13| (and references therein). 

In present work the necessary conditions for Fermion representability 



closely related to the Q condition of Coleman ||TT|, 12] are obtained for the 



most general case. Our technique is essentially finite-dimensional and ex- 
tensively use set-theoretical manipulations, enumerative combinatorics, and 
multilinear algebra. For these reasons we found it convenient to change 
notations that seems to be considered as traditional in literature on the rep- 
resentability problem. Namely, capital Roman letters are normally used to 
denote subsets of the spin-orbital index set N and/or the orbital index set 
M. Number of elements (indices) in index sets are usually denoted by small 
Roman letters, e.g. \N\ = n, \M\ = m, etc. With such notations it seems 
consistent to use small Roman lettes for the number of electrons and for the 
current density matrix order that are just the numbers of elements in the rel- 
evant spin-orbital index sets. Throughout this paper n is the number of spin 
orbitals, m is the number of orbitals, p stands for the number of electrons, 
and q is the density matrix order. 

The following results described in present paper we consider as novel: 

Explicit expression for A(n,p,q) operator is derived and its connection 
with Coleman's Q operator in the cases q = 1 and q = 2 is established; 

It is proved that A(n,p,q) in invertable and explicit analytic expression 
for its inverse is obtained; 

Distance function for the convex body giving an exterior approximation 
of the set of representable density operators is obtained; 

Analytic expression for the right inverse of the contraction operator as- 
sociated with the expansion operator is derived; 

Description of the exterior approximation for the set of all representable 
matrices is obtained in orbital basis for many electron systems with fixed 
total spin projection; 

Two algorithms for direct determination of approximate 2-density matrix 
are developed. 
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2. Contraction and Expansion Operators 



Let T n ,\ be one-electron Fock space spanned by an ortho normal set (ipi) 
of n molecular spin-orhitals. Electronic Fock space is defined as 

J~n — (J) Fn,p (1) 

where 

p 

Fn,p = /\Fn,l (2) 
Fnfi = C (3) 

where C is the field of complex numbers. 

"Determinant" basis vectors of the Fock space are conveniently labelled 
by subsets of the spin-orbital index set N: for any R G N the corresponding 
basis determinant will be denoted by \R). 

Creation-annihilation operators associated with spin-orbital index % are 
defined by the following relations 

a\\R) = (l-6 i>R )(-iy\RUi) (4a) 

a t \R) = M-l) W> (46) 

where 

x r l, if ieR ,s 

may be considered as a possible generalization of the Kronecker S symbol, 
and 

e=|{l,2,...,i-l}ni*| (6) 

is the sign counter. 

Let us introduce step-down and step-up (super)operators acting on the 
operator space T n ® T^- 

Cij : z — > diZOj (7a) 
Uij : z — > a\za,j (76) 
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where z is an arbitrary operator over T n - If z G T n ,r ®?ns then Cy-(z) G 
The operator space T n ®T^ may be equipped with the trace inner product 

(z\t) = Tr(zH) (8) 

It is easy to see that the step-down and step-up operators are mutually 
conjugate with respect to this inner product. 

Now we can define the contraction and expansion as 

c = J2cu (9a) 

i 

U = J2 U H ( 9& ) 

i 

Their restriction to the p-electron section of the operator space is iden- 
tical (up to the nonessential combinatorial prefactor) to the commonly used 
contraction and expansion. 

Let z p G T n . p ® J-'np be some p-electron operator. It is expanded via basis 
"determinant generators" as 

(p) 

zp = J2z RS \R)(s\ (io) 

R,S 

where the upper summation index (p) indicates that the sum is taken 
over all p-element subsets of the spin-orbital index set N. It can be shown 
H that 

( P ) (fe) 

c k (zp) = k\^Z R s £ (-lf RAS ^\R\K)(S\K\ (11) 

R,S KcRnS 

The definition and properties of the set- theoretical operation A are col- 
lected in Appendix A. 

Let us consider the electronic Hamiltonian of the form 

H = X^'WiX^ + o 12 (ij\ — \kl)a\a]aia k (12) 

i,j 1 i,j,k,l Tl2 
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which is obviously a particle number concerving operator: HT n ^ C 
for any p — 0, 1, 2, . . . The energy functional corresponding to some p-electron 
system is defined as 

E(z p ) = Tr(Hz p ) (13) 

and is a linear mapping from T n # <S> T* n v to the field of complex num- 
bers. In the particular case of a pure p-electron state z p = \^/)(^/\ Eq.(13) 
reduces to the standard average energy expression. Using specific form of 
the electronic Hamiltonian, it is possible to contract the energy domain and 
redefine energy in terms of g-electron operators (q>2). Introducing reduced 
Hamiltonian 

H p ^ q = -J2(i\ h \j)4 a j + o 2 (#1 — \kl)a\a)aia k (14) 

we can rewrite the energy expression (13) in the form 




The unknown operator on the right-hand side of this equation is q-electron 
one with the unit trace. Therefore, if we were interested in the energy func- 
tionals on jF n p <g> J 7 * we could easily redefine them on T nA ® T^ q with 
2 < q < p. Unfortunately, the actual situation is much more complicated. 
Indeed, physically relevant are the so-called density operators that are dis- 
tinguished from general linear operators by the following three properties: 

(i) Positive semidefinitness; 

(ii) Hermiteancy; 

(iii) Unit trace. 

The set of all p-electron density operators will be denoted by £„ iP . The 
set £n tP may be characterized both analytically and parametrically. Indeed, 
axioms (i)-(iii) give us analytic description of £ njP in terms of linear equalities 
and inequalities. On the other hand, it is not difficult to recognize that £„ iP 
is a convex hull of the so-called pure states: 

(;) 

t p e£ niP ^t p = ^2X i \^ i )(^ i \ (16) 
i=i 
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where 




(17) 

E^ = 1 (18) 



i=i 

and e T ntV . 

It can be easily verified also that 



^jC p q £ n , P C £ n ,g (19) 



The set 



Wn, P ,q = ^- q (£n, P ) (20) 

is a proper subset of £ nA and is called the set of all p-representable density 
operators of order q. Unfortunately, in general case the structure of these sets 
seems to be so complicated that numerous attempts to find their constructive 
description have not lead to practically useful results. 

Taking some inner point of W n)Pjg , we can try to construct functions that 
determine the distance from a fixed inner point of W n ,p,q to its border dW njPtq 
along some chosen direction. There exists natural central inner point of the 
convex sets £ n , q and W n ,p,q (the normalized identity operator): 



1 (i) 



w q 



jn\ E WW (21) 

where N is the set of spin-orbital indices, \N\ = n. With respect to this 
central point the energy functional may be rewritten as 



E(tq,e(t q )) = ^Tr(H p ^{w q + e(t q )[t q - w q ]}) (22) 

where e(t q ) determines the distance from w q to the border of W n , P ,q along 
the direction t q — w q . A technique of distance function e(t q ) values calculation 
involving minimization of a certain auxiliary function can be developed. In 
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our opinion, however, such an approach may be even more complicated than 
the commonly used CI one. Anyway, in present paper we will not try to 
follow this route. 
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3. Exterior Approximation for the Convex Set 
of Representable Density Operators. 



First of all let us analyze one of two cases where the ensemble repre- 
sentability problem admits constructive solution ||. To this end it is con- 
venient to introduce a new basis in the operator space that differs from the 



basis of the determinant generators only by the phase transformation [ 14 1 : 

e[ J = (-l)K /UJ ) nAi l|/UL)(JUL| (23) 
where I D J — 0. From this definition and Eq.(ll) it readily follows that 

(?) 

tf^4 J = (P (24) 

KcL 

(?) 

u p - q e^ = (p-q)\ £ 4 J (25) 

Over the finite-dimensional electronic Fock space there exists an involu- 
tion defined as 

l\R) = \N\R) (26) 

Its restriction to the p-electron section of the Fock space is an isomorphic 
mapping from T UiV to T n ^ n - V and 

Vn- P = idjr nn _ p (27a) 

I n _ p I p = id Tnp (276) 

where id^ n n _ p and idp are the indentity operators over T n ,n-p and J r n ,p, 
correspondingly With each vector \^f) G T n $ a semilinear functional (^/\ : 
^n,p — > C is associated and 

(*|ln_p = (28) 

where {^\l n - p is a composition of mappings T n , n -p — Fn,p ~^ C. 
For each pair (n — p, p) let us put 

\i—p,p • Zp * ^-pZp^n—p (29) 
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where z p is an arbitrary p-electron operator. It is immediately clear that 
I n _ PiP is an isomorphic mapping of <g> !F* V onto T n ,n-v ® F* n_ p . Simple 
manipulations lead to the conclusion that 

\i—p,p(.£n,p) £n,n—p (30) 

The effect of l n -p, P on the basis operators (23) is given by 

i n - P M J ) = (-in j 4 7 (31) 

where a u = \(I U J) n A N \ {IUJ) \, and K = N\(I U JUL). Note that 
the sign prefactor on the right-hand side of Eq.(31) depends on the index set 
I U J but not on the number of elements in it. 

Let us introduce certain symmetric combinations of g-electron and p- 
electron basis operators: 

™&W = 7^X?# (32) 

q-s) K ^ L 

(p-s) 



<U K ) = T^A E 4 J (33) 



n—q— 
n—p—s 



where I n J = 0, |/| = | J\ = s, K, L C N\(I U J), \L\ =p - s, and 
\K\ = q — s. 

Using definitions (32)and (33), we can rewrite Eqs.(24) and (25) in the 
form 



<t^, e u = ^k <lq{L) (34) 

Let us suppose that p + q = n. In this case the contraction operator is an 
isomorphic mapping || [TJj] which means that the set of vectors (32) is a basis 



ji 



set of the q-electron operator space. We can therefore define nondegenerate 
linear operator in the following manner: 
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A(n,p, q)w I p J lq (L) = (-ir iJ if^ e J N \ (IUJuL) (36) 



p 

q-s 



Comparing Eqs.(31),(35),and(36) immediately leads to the conclusion 
that for p + q = n 

A{n^q)^(f-^l n ^ p (37) 



and, consequently, 

A(n,p,q)W n , p , q = £ niq ^ 



(p + q = n) 

It seems pertinent to note that operator satisfying Eq.(38) is not unique. 
We can as well take any composition of A(n,p, q) with operator leaving the 
set £ n ^ q invariant. For example, the transposition 

r:\R)(S\^\S)(R\ (39) 

possesses this property and we can use rA(n,p, q) instead of A(n,p, q). 
In contrast to the operator A(n,p, q) the operator rA(n,p, q) is conveniently 
block-diagonal: 

rA(n,p,q) = ($tA ij (n,p, q) (40) 
i, J 

where A IJ (n,p, q) is the restriction of A(n,p, q) on the subspace of p-electron 
operator space spanned by the basis operators (23) with fixed /, J. 

After complicated combinatorial manipulations (closely related to the fa- 
mous inclusion-exclusion principle) the explicit matrix representation of the 
operator rA(n,p,q) can be obtained (see Appendix B) 



(P) (q-s) (p+\KnK'\-q-l\ 

T ^ 4 )ek' = (-ir' 7F 4«U e (-Di W ' ex 

[q-s) { q ) K'CN\(IUJ) \\KnK'\) 

(41) 

By an abuse of notation, we will use the symbol A(n,p, q) for the operator 
defined by Eq.(41). 
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Thus, in the case p + q = n the convex set ~W n ,p,q can be explicitly char- 
acterized in several equivalent ways: 

as the pre-image of E nA with respect to A(n,p, q) 

W n , Ptq = A- 1 (n,p,q)S ntq ; (42a) 

as the set of solutions of the system 



({$ q \A(n,p,q)t q \$ q )>0,$ q eF n , q 
\Tr(t q ) = l; 

as a convex body with the distance function 



(426) 



e(t q ) = — (42c) 

mm ($ q \A(n, P ,q)t q \$ q ) 

1*911 = 1 

where $ 9 G T n , q and Tr{t q ) = 1. 

Let us try to generalize the approach described to handle the case p + q < 
n. The first idea coming to mind is to use the operator defined by Eq.(41) in 
general case. This operator is nondegenerate for n > p > q, p + q < n and 
explicit analytic expression for its inverse can be obtained (see Appendix C). 
It can also be shown that 

A(n,p,q)W niPiq C £ n , q (43) 

(see Appendix D). From this inclusion and positive semidefiniteness of 
the reduced density operators it readily follows that the convex set 

V n , p , q = A-\n, p, q)S n<q n £ n , q (44) 

is an exterior approximation for W n , P , q - Simple combinatorial manipulations 
lead to the conclusion that 

A(n,p, q)w q = w q 

that is V n , p , q , Wn, P ,q, and £ n ^ q share the same central point. 

Thus, we have constructed the compact convex set that may serve as a 
certain exterior approximation for the set Y\? n ,p,q and coincide with W njPig in 
the particular case p + q = n. It is not difficult to demonstrate, however, 



11 



that in the case q — 1 this set differs from W n ,p,i- The reason is the unitary 
noninvariance of A(n,p,q). 

It is well-known that the contraction operator is invariant with respect 
to the induced unitary transformations of the type 

k k 

p k (u) = (/\u) <g> (A^), u G U(F nA ) (45) 

that is 

cp k {u)t k = p k _i(u)ct k (46) 

for any t k G T r ^ k ® J~n,k- Taking into account obvious unitary invariance 
of the set £ njk , we can conclude that 

p q (u)W n , p , q = W n , Pig 

for any u G U(J 7 nj i). 

Direct calculations show that to restore the desired unitary invariance 
of V n ,p,i h is sufficient to modify the phase prefactors in definition (41) (see 
Appendix E). 

Let us introduce the phase transformation 

v : \R)(S\ -> (-iy RAS)nANl \R}(S\ 

It leaves invariant the sets S n>q and, consequently, we can replace the op- 
erator A(n,p,q) by uA(n,p,q). Matrix representation of the last operator 
differs from Eq.(41) only by the sign prefactor that should be taken equal 
to (— l)K /UJ ) nA (J"u./)l i Hereafter, only the operator uA(n,p,q) will be under 
consideration and, by an abuse of notation, we will use for it the same symbol 
A(n,p,q). 

In the particular case q = 1 we have 

A(n,p,l)e» = ±l^L e f ( 47 a) 
i=i n ~P 

A(n,p,l)e$ = V —eX (476) 

n — p 

Using the unitary invariance of A(n,p, 1) we can state that t\ G V niPi i if 
and only if uH\U G V n ,p,i where u is the unitary transformation diagonalizing 
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density operator ti.But for diagonal t%, with the aid if Eq.(47a), we can 
obtain the equivalence 

t 1 = Y t ^eV n , p , 1 ^l [ ° T f { ^- P (48) 

that immediately implies V n , Pj i = W n , P ,i (consequence of the diagonal ver- 
sion of the famous Coleman's theorem Q). 

For the cases q=l and q=2 it is easy to recast operator A(n,p, q) in more 
invariant form. We have 

~ P j A{n, p, l)ti = - Q h + Tr(t x )id Fnil (49) 

and 



2 P ^A(n,p,2)t 2 = ^V2-p(*iAi^ nil +i^ nil Ati)+Tr(* 2 )A^n,i (50) 

where idj? nl is the identity operator over jF n l and t\ = \ct 2 is 2 — ^ 
1 contracion of t 2 . The operator (50) is similar but not identical to the 
Coleman's operator Q |JTT], TJJ. Operator A(n,p, 2) differs from Q in a 



rather delicate manner. Namely, 2 p ^jA(n,p, 2)t 2 includes symmetric term 
p{ti^\id Fnl +id Fnl Mi) whereas Qt 2 contains the term 2p(t\f\id Fn J which, in 
our opinion, is not correctly defined. The situation here is completely analo- 
gous to that in the angular momentum theory. For example, for two-electron 
systems one-electron operator of momentum projecton j z has no sense in it- 
self. Only the total operator J z = j z Aid Fnl + id Fnl A j z is of mathematical 
and physical meaning and it can not be replaced by 2(j z A id Fn l ). 

The sets V n , P , q may be explicitely described in several equivalent ways: 
As the intersection of £ n>q with its pre-image with respect to A(n,p, q)(see 
Eq.(44)); 

As the set of solutions of the system 

(%\t q \%) > 0, 

{$ q \A(n,p,q)t q \$ q ) >0 (51a) 



I Tr(t q 



where $ g e T n ^ 
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As a convex body with the distance function 



e(tg) = min{ ^ , ^ } 

1 - (^mmj%\t q \%) 1 - ($mmJ$ q \A(n,p,q)t q \$ q ) 

(516) 

where <E> g G J- n>q and Tr{t q ) = 1. 

From the well-known properties of the distance functions (see, e.g., ||15|| ) 
it follows that t q e dV n ^ q if and only if e(t q ) = 1. But from Eq.(51b) it is 
easy to see that e(t q ) = 1 if and only if either t q or A(n,p,q)t q has at least 
one zero eigenvalue. In the particular case of pure p-electron determinant 
state |-R)(-R| the corresponding g-density matrix 

i (?) 

f^) = 7ST E l*)<*1 (52) 







has — zero eigenvalues and, being obviously representable, belongs 
to 9V n , Pi g n dW n ^ q . 

Direct but somewhat tedious combinatorial calculations lead to the im- 
portant commutation relation of the operator A(n,p, q) with the contraction 
operator (see Appendix F): 

A(n,p, q — l)c = cA(n,p, q) (53) 

Now we can collect the most important properties of the convex set V niPiq : 

(1) It gives an exterior approximation for the set W n , P ,q of representable 
density operators; 

(2) From Eq.(53) it follows that 

V n , P ,q = yVn, P ,q for q = n-p and g = l; (54a) 

^cV n , p ,g C V n , P ,q-i for any q <p (546) 

(3) The set V n ,p,q is invariant with respect to transformations induced by 
unitary 1-electron ones (see Eq.(45)); 

(4) All Hartree-Fock densities belong to the border ofV n)Piq ] 
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(5) Any representable q-density matrix t q such that either t q or its image 
with respect to A(n,p,q) possess zero eigenvalue belongs to the border of 

Density matrix of order q associated with pure p-electron state practically 
always possesses property (5). Indeed, if there exits one-electron basis with 
at least one molecular spin-orbital having occupancy (virtual) or 1 (inac- 
tive) then block-diagonal matrix t q © A(n,p, q)t q necessarily has zero or very 
small eigenvalue. Practice of ab initio calculations shows that even in bases 
of moderate size among natural spin-orb itals there present inactive and/or 
virtual ones. We may therefore hope that the direct energy optimization on 
the set V n ,p,q will lead to reasonable results. 

For readers who prefer numerical arguments instead of abstract mathe- 
matical ones, in Section 5 the results of relevant calculations testing property 
(5) of the set V niPi2 are presented. 
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4. Sections of the Contraction Operator 

The contraction operator p jC p ~ q is linear surjective mapping from jF n p <g> 
J-* on T n , q ®J~n q - By definition, global section of the contraction operator 
(or its right inverse) is linear injective mapping ir q ^ p from T nA ® T* n q to 

Fn,p ® Fn,p SUCn ^at 

-C^ % Tp = %d Tn ^ Tlq (55) 

Here we describe a certain general scheme of construction of the con- 
traction operator sections and illustrate this scheme on example of section 
associated with the expansion operator.. 

First of all we should find an injective linear mapping 7^ of the form 

W40= Y; 4 J f"(L,K) (56) 

where is a certain set of (p — s)-element subsets from N\(I U J). 
Injectiveness of the mapping (56) means that vectors 7 ? t P (e^) are lineary 
independent. Then we should calculate 

4^7^)= ( E &G% K (57) 

K'CN\(IUJ) 

where we introduced (in general non-symmetric) matrix G(n,p,q) with 
matrix elements 

G^ = 4 ( E (58) 

If this matrix is invertable then the desired section of the contraction 
operator may be written as 

n q ip = l q ipG~ 1 (n,p,q) (59) 

As an example, let us consider the section associated with the expansion 
operator (see Eq.(35)): ^ q ^ p = j^z^u p ~ q . In this case in Eq.(56) we have 
= {L C N\(I UJ):LDK}, and f IJ (L,K) = -J^-. First of all we 



\p—ij 



16 



should ascertain that the expansion operator is injective. To this end it is 
sufficient to show that the vectors (33) are lineary independent. Then we 
should construct matrix G(n,p,q) and try to invert it. All relevant technical 
details are collected in Appendix G. Here we just write down the final explicit 
expression for the section under discussion: 



(p) (p-s) (p-s-\KnL\-l\ 

k)=(-^- s 7^u e 4 j (-i) ixnLi V^V ( 6 °) 

LcN\(IUJ) \q-s-\KnL\) 



p-q 



For the case p + q = n there exists the unique global section of the 
contraction operator equal to its inverse (see Eqs.(29) and (41)): n q ^ p = 
R p ,n-pA(n,p,q). We can try to generalize this last section to treat the case 
p + q > n. One of the ways is to define vectors 



(P) (p-s) fp-s-\KnL\-l\ 

= (-ir'j^Lx s e"(- 1 ) ,LnKI T*--\ (61) 

\qs){ q ) LCN\(IUJ) \\KnL\) 

and consider the mapping 7 ? | p (e^) = . 

Besides global sections satisfying Eq.(55) of interest are also sections 
closely related to the concrete density matrix under consideration and de- 
fined on a certain subspace of T n ,q ® They can be used to restore the 
part of p-density matrix that survives under contraction to the level q. Any 
progress in this direction would be of great importance. 
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5. Density Operators in Orbital Representation 



Let (ifi) be an orthonormal set of m orbitals. Following Handy |[16 



we 



identify p-electron determinants generated by these MOs, with pairs of index 
sets (strings): 

\R a , Rp) = a\ ia . . . a\ vaCt a) iP . . . a) ppP \Q) (62) 

where R a = 1 < i\ < . . . < i Pa < m, Rp = 1 < j 1 < . . . < j p < m, 
Pa + Pp = P, and |0) is the vacuum vector. 

Split basis operators analogous to ones defined by Eq.(23) are 

where I a fl J a = Ip fl Jp = 0, 

ei = \I a UJ a \ x \Lp\ 2 + \(I a UJ a )nA La \ + \(IpUJp)nA Lg \, 



2 



0, if \Lp\ = (mod 2) 

1, if \Lp\ = 1 (mod 2) 



and 



( i^n-qA I a,Ja){Il3,Jf3) _ 1 \ ^ \ ^ V~» {la ,Ja ) (1/3 ,Jp) 



ka,*0 K a CL a KpCLg 

(ka+k^—q — s) 



. ro (J-. J »)('«"«)(L„, L(J ) (64) 
The operator A(2m,p, q) acts on the basis vectors (63) as 



A(2 m , P ,,)e&S^ J » , = (-l)« 7 -J—V £ E 



1 2m— p 



q J k a ,kp K' a CM\(I a UJ a ) 

{fea+fc^ — q — s) 



{k } /p+|K a n^|+|^n^|-g-i\ 

/ 1 x|g a n^|+|^ fl n^| V l^n^|+|^n^| ^ (I a ,J a )(l(,,Jp) (Q5) 
2-~t ' ( as \ K,^) \ ' 



K' CM\(I UJp) 



{\K a nK'J+\KpnK' \) 
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where M is the orbital index set, s = \I a \ + \Ip\ = \ J a \ + \ Jp\, and 

e 2 = \Ia U J a \ X \Ip U Jp\ 2 + \(I a U J Q ) H A (JaUJa) | + | (/^ U Jp) fl A (J<jUj/3 )| 

Note that the sign counters ei and e 2 in these formulas are determined by 
the initial spin orbital index set ordering and its current form corresponds to 
the split determinant representation (62), where a indices always go first. 

Let us consider p-electron states with a given value of the total spin 
projection M$- The set of all such states is generated by the determinants 
(62) with fixed |i? Q | = p a , \Rp\ = pp, since, by assumption, p a + pp = p, and 
Pa—pp = 2M5. Contraction of an arbitrary determinant generator p— 2 times 
may lead to nonzero result only if \R a C\ Rp \ + \ S a n Sp\ > p — 2 which means 
that on the 2-electron level there may appear only three types of generators: 

|rir 2 ,0)(sis 2 ,0| 
|0,r 1 r 2 )(0,s 1 s 2 | 
\n,r 2 )(s 1 , s 2 \ 

As a result, of interest are block-diagonal 2-electron density operators of 
the form 

t 2 = t% + 4 + tf (66) 



where 



*2 = E ^ lS2 |nr 2 ,0>( SlS2 ,0| (66a) 

r 1 <r 2 
s 1 <s 2 

4= E C 2;siS2 |0,rir 2 )(0, SlS2 | (666) 

r l <r 2 
s 1 <s 2 

tf = E ^r 2;SllS2 ki,r 2 )(si,s 2 | (66c) 



ri,r 2 



Calculation of traces of 2-density matrix components corresponding to 
some p-electron wavefunction of fixed total spin projection gives 

UJ (67a) 
(<r = a,P) 
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and 

Tr{tf) = ™* (676) 

Eqs.(67a)-(67b) may be considered as additional restrictions on density 
matrix components in the course of energy optimization. 

Using either general definition (65) or invariant operator (50) one can 
show that the action of the operator A(n,p, 2) on arbitrary 2-density matrix 
t 2 of block-diagonal form (66) results in matrix of the same block-diagonal 
structure 

d 2 = A{n, p, 2)h = d% + d£ + df (68) 

and trace restrictions (67a)- (67b) imply the following restrictions on the 
traces of spin-components of operator d 2 . 



( 2 T P ) 



[a = a, /3) 

and 



(69a) 



Tr(df) = im -^p-^ (6%) 



Reduced Hamiltonian in the orbital representation is 
1 A , „ 1 



p^2 — — 
P 

where 



+ 2 E (immjEkt-SjkEu] (70) 

i,j=l i,j,k,l=l 



h ij = J Viir^hipjir^dn (71) 

R 3 

(ij\kl) = f (p^rjtpjin)— (p* k (r 2 )(pi(r 2 )dridr2 (72) 

R 3 XR 3 

and 

E ij = E ij + E ij = a ia a ja + alp a 3f3 ( 73 ) 
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are the unitary group generators. 

Standard but somewhat tedious algebraic manipulations lead to the fol- 
lowing general energy expression 



E(h) 



P 



m m 

EE 

i k=i 



(t a . _)_ jj 3 . + t al3 ■ + 1 



.a/3 



)]hij+ 



+E(*S;H + *5;H)[(* i iy)-(^l K )]+ E 

i,j,k,l=l 



k<l 



(74) 



where 



L ij-,ki — 



v ij;kl 
L ij,lk 
L ji,kl 
L ji;lk 



(75) 



if i < j and k < I 
if i < j and k > I 
if i > j and k < I 
if % > j and k > I 
if i — j or k — I 
(a = a, (3) 

The property (5) of the exterior approximation for W2 m , p ,2 given by the 
convex set V2 m , P ,2 (see Sec. 3) may be tested numerically in the following 
rather obvious manner. With the aid of CI or CASSCF method one can ob- 
tained p-electron wavefunction expansion over determinant basis set and then 
contract the pure p-electron state to get spin components of the 2-density 
operator. Calculation of the lowest eigenvalues of six 2-electron operators 
t% , if, W i $2 ) d-2 , and with subsequent use of Eq.(51b) to get £(£2) makes 
the problem of energy E{ti,e(t^)) (see Eq.(22)) evaluation at the border of 
V2m,p,2 very simple. It is pertinent to note that we can take either the normal- 
ized identity matrix as a central point of the set V2 m , p ,2 or its "spin-adapted 
analogue" : 



w 2 



<P<*\ Pf3\ 



T^EI^^^^I + T^EI^^^^I + ^EKi)^^ 

(76) 

For energy evaluation at the border points of V2 m ,p,2 the choice of central 
point is of no consequence. At the same time for optimization purpose central 
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point denned by Eq.(76) should be used to guarantee that at each step of 
energy minimization current 2-density matrix has block-diagonal form (66). 
Note also that selection of the central point (76) requires slight modification 
of the expression for the distance function e(t 2 ) of the convex body V2m, P ,2- 
Algorithms for contracting CI expansions to get the corresponding 2- 
density operator components t^t^a-iad t% and subsequent transformation 
of these components with the aid of the operator A(2m,p, 2) were devel- 
oped. Using these algorithms in parallel with GAMESS program set [17 , we 



performed calculations of CI wavefunctions for ground and excited states of 
small atomic and molecular systems to get for each pure state the distance 
function £(£2) value to estimate its proximity to the unit and to compare the 
energies at the border of V2 m ,p,2 with the CI energies. The detailed scheme 
of the algorithm used can be described as follows. 

(1) Calculate CI wavefunction \1/ in determinant basis set; 

(2) Contract pure state (\&| to get spin components of 2-density matrix 

(3) Calculate spin components of the matrix A(2m,p, 2)t2(^); 

(4) Evaluate lowest eigenvalues of six 2-electron matrices t%(ty), d?,^) 
where a = a,/3,a/3, and calculate the distance function value efo^)); 

(5) Calculate energy of the central point used; 

(6) Calculate energy value at the point w 2 + £(t 2 (^))[t2 — w 2 ] G <9V2 m , P ,2- 
Note that the final distance function and energy values are of no primary 

interest because they should be close to the unit and the CI energy value, 
correspondingly, due to the property (5) of the convex set V2 m , P ,2- From 
Eq.(51b) it follows that the actual distance to the border of V2 m , P ,2 is deter- 
mined as the minimum of distances from the central point to the border of 
^2m,2 and the border of A~ 1 (2m,p,2)£2m,2- Of certain interest is which of 
these two distances determines the final one for concrete pure CI state. 

The results of trial calculations are listed in Tables 1 and 2. In atomic 
calculations cc_ pvDZ basis set of Dunning |18[] was employed. For lithium, 
berrilium, and boron FCI calculations were carried out, whereas for carbon 
and nitrogen Is AO, and for oxygen Is, 2s AOs were excluded from the active 
spaces to keep the sizes of CI expansions reasonable for running GAMESS 
on PC with 166 MHz Intel processor. In the case of molecules cc_pvDZ 
basis and FCI were used for calculation of LiH, and 6-31G Gaussian basis set 
U, ||] with frozen Is AO for CH 2 , and Is, 2s AOs for NH 2 ,H 2 0, and NH 3 
was employed. 
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As seen from Tables 1 and 2, the exterior approximation for pure CI states 
given by the convex set V2m, P ,2 is very good as it was expected. It is to be 
noted that the contraction procedure may lead to roundoff errors in matrix 
elements of 2-density matrices, especially when CI expansion is large. As a 
result, diagonalization of 2-density matrix may give small negative values for 
its lowest eigenvalues in spite of the fact that the matrix t 2 (^) for any ^ is 
manifestly nonnegative. For this reason in Tables 1 and 2 the absolute values 
of e(t 2 ) deviation from the unit are given. Of interest is also the fact that 
for atoms and high symmetry linear molecules the distance from the central 
point to the border of V2 m , p ,2 is determined by the lowest eigenvalue of the 
relevant density operator whereas for low symmetry molecules this distance 
is determined by the lowest eigenvalue of the operator A(2m,p, 2)t 2 - 
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6. Algorithms for Direct Determination of 2-Density Matrix 



The energy expression (74) may be essentially simplified by turning to 
a new basis set in the 2-electron section of the Fock space. Indeed, let us 
consider the basis set of eigenvectors of three matrices 

ui ij;kl 

Km = (™) 

UL ij;kl 

p _ dE{t 2 ) 

U i,j;k,l ~ r. ap \'< C ) 
Ul i,j;k,l 

that are just (up to fixed prefactor) spin blocks of the reduced Hamil- 
tonian H p ^ 2 - In this "energy" basis arbitrary 2-density operator is of the 
form 



h = E Ks\^)m + E W)(<&f| + E K!\*f)(*?\ (78) 

r,s—l r,s=l r,s=l 

and the electronic energy expression may be written 

\2 ) \2 ) m 2 

E{h) = E + E + E (79) 

r=l r=l r=l 

Here {e", |$")}, {e^, |$f )}, and {e"^, l^^ 3 )} are the eigenvalues and eigen- 
vectors of the matrices (77a), (77b), and (77c), correspondingly 

Since in the energy representation the electronic energy is a linear func- 
tional involving only the diagonal elements of 2-density matrix it seems rea- 
sonable to try the well-known Box method [21] for energy optimization. This 
method starts with generating at least k = ( 2 ™) + 1 affine independent di- 
agonal matrices A^ 1 ), \( 2 \ . . . , \( K > with nonegative entries such that for each 
diagonal matrix A^ from this set the following Approximate Represent ability 
condition (yli?-condition) is satisfied: 



24 



(AR) There exists (at least one) symmetric matrix with zero diagonal 
entries such that 

U\®W + UA®W e V 2m , p , 2 (80) 

where U is the unitary transformation from the 2-electron determinant basis 
to the energy basis. 

It is pertinent to note here that the convex set V 2m)P)2 being invariant with 
respect to the orbital unitary transformations is actually changed under 2- 
electron transformations. 

Then we should order the initial vertices in the energy increasing order: 

E(\ (1) ) <E(\W) < ... <E(\W) (81) 
and calculate the weight center of the first k — 1 vertices 

A~D (0 (82) 
K ~ 1 i=i 

This weight center obviously satisfies the Ai?-condition and is used in the 
Box method for reflection of the "worst" vertex A^: 

A( e ) = A + e(A- A (K) ) (83) 

where e is a positive reflection coefficient. It is clear that E(X) < E(X^) 
and that the energy value may only decrease when moving along the vector 
A — \( K \ The optimal value e* of the reflection coefficient should satisfy the 
following condition: 

(AR) There exists (at least one) symmetric matrix A with zero diagonal 
entries such that 

U\(e*)U^ + UAW E dV 2m , p , 2 (84) 

where U is the unitary transformation from the 2-electron determinant basis 
to the energy basis. 

The next step of the Box algorithm consists in replacing the "worst" 
vertex by the calculated optimal vertex A(e*) with subsequent reordering 
the new set in the energy increasing order. Then the new weight center should 
be constructed and new optimal vertex should be determined. This process 
is repeated till the energy value stabilization and the final weight center is 
taken as the desired solution. 
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The crucial point of the algorithm described is the calculation of the 
optimal vertex satisfying the Ai?-condition. Let us suppose that the current 
2-density matrix is of the form 

t 2 =t x + UAU^ (85) 

where diag(UH x U) = A, and A is some symmetric matrix with zero diag- 
onal entries. The distance of matrix (85) from some symmetric nonnegative 
2-density matrix X 2 (in the standard trace metric) is given by the function 

/j 1} (A, X) = Tr(t x + UArf - X 2 )\t x + UArf - X 2 ) (86a) 
whereas the distance of this matrix from some matrix A^ 1 (Y 2 ) is 



/f )(A, Y) = Tr(t x + UAU* - A- l {Y 2 ))\t x + UAU ] - A^iY 2 )) (866) 

Here X and Y are symmetric matrices of parameters that have the same 
block-daigonal structure as 2-density matrix (66) and satisfy the restrictions 
on block trace values given be Eqs.(67a)-(67b) and Eqs.(69a)-(69b), corre- 
spondingly. 

The function 

fx (A, X, Y) = ( A, X) + /f (A, Y) (87) 

can be considered kind of characteristic function of the set V2 m v 2 

of all approximately representable 2-density matrices. Indeed, if t 2 G V2m, P ,2 
with diag(UH2U) = A then 

f x (0,tl[A(2m,p,2)t 2 ]^)=0 

On the other hand, if for fixed t x there exist matrices A, X, and Y, such 
that 



f x (A,X,Y)=0 

then, due to the nondegeneracy of the scalar product used, t 2 = t x + 
UAW = X 2 = A- 1 (2m,p,2)(Y 2 ) belongs to V 2m , P ,2- 

Now we can formulate the (AR) condition in the form admitting numer- 
ical verification: 
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(AR)' If for a fixed diagonal 2-density matrix A there exist symmetric 
matrix A with zero diagonal entries and symmetric matrices of parameters 
X and Y such that 

f x (A,X,Y)=0 

then 

UXU^ + UAU^ e V 2m , p , 2 

To calculate the optimal vertex A(e*) it is necessary to analyze the non- 
negative scalar function 

f(e)= mm Y f m (A,X,Y) (88) 

that is strictly positive for e > e* and accepts zero values in the interval 
[0,e*]. 

Note that the initial set of vertices in the Box method should be carefully 
selected to ensure their affine independence because in the opposite case, 
due to linear character of this method, the search of optimal vertex will 
be performed only within a face of the polyhedron spanned by the initial 
vertices. 

The method described being conceptually very simple is not easy for 
practical implementation. Its expected main drawbacks are: 

(1) It requires often transformations from 2-electron determinant basis to 
the energy one and back; 

(2) Optimization algorithm include both differentiable and non- differen- 
tiable parts; 

(3) Minimization of function (87) with respect to X and Y should be 
performed with very high accuracy (about 1CT 20 in gradient norm); 

(4) Function (87) has excessive number of parameters that results in de- 
generate Hessians and low convergence of minimization process. 

There exists another approach based on the energy expression modifica- 
tion. Since operator A(n,p,q) is symmetric and invertable (see Sec. 2 and 
Appendix C), we can write down the following equality: 

Tr[H p _> 2 t 2 ] = Tr[(A- 1 (2m,p, 2)H p ^ 2 )A(2m,p, 2)t 2 ] 

Introducing two symmetric matrices of independent variables we can 
rewrite the energy expression (22) in the form 
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(p) 

E(X,Y) = ^-(Tr[H p ^ 2 X 2 }+Tr[(A- 1 (2m,p,2)H p ^)Y 2 } ) (89) 



2 

The corresponding optimization problem 

' mm E(X,Y) 



X.Y 

-2 



A(2m,p,2)X 2 = Y 2 m 
Tr{X 2 ) = 1 
{ TriY 2 ) = 1 

may be solved by constrained optimization technique (see, e.g. p2| ). 
This approach, in contrast to the first one, does not require 2-electron uni- 
tary transformations and both the energy expression (89) and the restrictions 
on variables are differentiable. Its main drawback is shared by all methods 
based on algorithms involving Lagrange multipliers and penalty functions: 
in theory infinite number of intermediate unconstrained optimizations is re- 
quired to reach the solution. It is to be noted as well that the number of 
variables in Eq.(90) can be essentially reduced if we take into account that 
matrices X and Y should be of block- diagonal structure for states with fixed 
total spin projection. If block-diagonal structure of X and Y is explicitly 
accounted, it is reasonable to replace trace resrictions for X 2 and Y 2 by six 
trace restrictions for the corresponding spin blocks (see Eqs.(67a)-(67b) and 
(69a)-(69b)). 

6. Conclusion 

The exterior approximation for the set of all representable density oper- 
ators of arbitrary order described in this work is expected to be useful for 
direct 2-density matrix determinantion. However, only on the base of con- 
crete calculations it will be possible to estimate the practical importance of 
the approximation obtained. 

Two algorithms for such calculations are developed. In contrast to the 
standard approaches where the electronic energy domain turns out to be one 
of classic analytic manifolds (unit sphere, orthogonal group and its quotients, 
etc) that can be easily parametrized by, say, elements of relevant tangent 
spaces |p3fl , the convex set V2 m , P ,2 is of much more complicated nature. The 
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only thing we can try do at present is to cover V2m,p,2 using excessive set of 
parameters. Energy optimization on V2m,p,2 may be a very complicated com- 
putational problem but if efficient optimization scheme is developed, there 
may be opened a way to FCI quality calculation of fairly extensive molecular 
systems. 
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Appendix A. 

Let N = {1, 2, . . . , n} be the spin-orbital index set. On the set V(N) of 
all subsets of N let us consider the operation 

RAS = (RU S)\{RnS) (Al) 

where R, S G V(N). This operation endows V(N) with Abelian group 
structure with empty set as its unit. Each element of this group is of order 
2 (RAR = 0). The mapping 

y? A : K -> A K (A2) 

where K C N and 



A*: = A fc6A :{l,2,...,A;} (A3) 

is a group homomorphism. Indeed, A0 = and (Ax)A(Al) = Akal- 
In particular, 

A W = {1,2,...,A;} (A4a) 

f{A; + 1,...,/} ifk<l 
\{Z + l,...,Jfe} ifk>l 

{2, 4, . . .} if n is even 
{1, 3, . . .} if n is odd 

In general case, for K = k\ < k 2 < . . . < k s 



A {k,l} = 

A N = 



A K = < 



■ z J 

U {hi-i + !,•••, hi} if s is even 



i=i 
[f] 



U {hi + 1, • • • , hi+i} if s is odd 

I i=0 



(A46) 
(A4c) 



(A5) 



Directly from the definition of operation A the following relations impor- 
tant for phase prefactors evaluation may be obtained 



\KDR\ + \KnS\ 

\KnA K \ = [ 



Kn(RAS)\ {mod 2) 
K\ + 1, 



(A6) 
(A7) 
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Appendix B. 

From Eq.(30) it is easy to get matrix representation of the operator 
A _1 (n, p, q) with respect to the basis {e^} for p + q = n. Confining ourselves 
to the diagonal case I = J = 0, we obtain 

f T^r, if K> C N\K 

A K ) K (n,p,q) = \ ©' (5.1) 
lO, iiK'<£N\K 

Let us suppose that 

A KK ,(n,p,q) = (-l)\ KnK y(p,q,\KnK>\) (B.2) 

A(n,p,q) matrix elements are determined by the equations: 

(?) 

J2A KK/ (n,p,q)A- K ) K „(n,p,q) = 

K< 

i (?) 

= 7^ E (-ij^/fegj^n^'l)^^ (5.3) 

y k>cn\k» 

After simple set-theoretical manipulations we come to the following com- 
binatorial equations 

s ( - i)r G : 3(9 /( ** r) =GK (B - 4) 

where t = \K\{K fl X") | and r = |if D if |. From this equation it is easy 
to get reccurently the expression for AKK'(n,p,q): 

(p+\K<lK'\-q-l\ 

A KK ,(n,p,q) = (-l)^ iip^ (5.5) 
Note that Eq.(B.4) can be rewritten in the form 

s-CHr:)^— (:)(';') ™ 

and proved by induction for p > q > 1 and £ = 0, 1, . . . ,p. 
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Appendix C. 



Let us reduce the problem of operator A(n, p, q) inversion to combinatorial 
equations. The matrix elements of the operator A(n,p, q) with respect to the 
operator basis (ej/) are given by Eq.(43). Since this matrix is block-diagonal, 
and all its blocks are similar in their structure, we can confine ourselves to 
the block with I = J = 0. 

Let us consider the system of linear equations 

(?) 

H A KK'( n ^P^) X K'K"(n,p,q) = 5\ K nK"\, q (CI) 

K'cN 

for determining the inverse matrix (block). Simple combinatorial argu- 
ments together with the additional assumption that Xk>k" = f( n ,P, <?; \K'n 
K"\) (that is X K 'K" depends not on subsets K' , K" but only on the number 
of elements in their intersection) allow us to rewrite the system (C.l) as the 
following system of combinatorial equations 



= ("~ P )^ (C2) 

where u = \K D K"\ = 0, 1, . . . , q, r\ = \K D K'\ = 0, 1, . . . , q, r 2 = 
\K' n K"\ = 0,1,..., q, and t = \K (~) K' (~) K"\ = 0,1,..., u. Thorough 
analysis of this system leads to the conclusion that 

/ n-p-q+r 2 -l\ 

f(n,p,q;r 2 ) = (-ir ± ^ > (C.3) 

and, as a result, the operator A(n,p,q) is invertible and the matrix ele- 
ments of the inverse matrix are 



(n-p\ (n-p~q+\KC\K'\-l\ 
V 5-* )\q) \\KnK'\) 
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where |7| = | J\ = s, I n J = 0, K, K' C iV\(J U J), |X| = \K'\ = q - s. 
The last equation implies, in particular, that 

= A(n,n-p,q) (C.5) 



33 



Appendix D. 

Lemma. 



A(n,p,q)W n , p , q C£ n , q (DA) 

Proof. It is sufficient to show that A(n,p, q)^c p ~ q \^} (^| G E nA for arbitrary 
^ G J-'n^p. Expanding \1> over p-electron determinant basis set 



(p) 

* = e c r \r) 

RcN 



and turning to the operator basis (e£ J ) in T U)V ® p we obtain 



^W(*I = EE I? (-i) l(/UJ)nAil ^uLC} UL iyl^(^) 

P- s=0 /,J LCN\(IUJ) [g) 

(D.2) 

where symmetric combinations Wpj q (L) of q-electron basis operators are given 
by Eq.(32). Simple set-theoretical and combinatirial manipulations lead to 
the equality 



\q-sj { q ) KCN\(IUJUL) 

Turning back to q-electron determinant generators and carefully handling 
the phase prefactors arising, we arrive at 

n \ 1 (n-p-q) (p) 

A(n,p,q)±<r*\*)(*\ = 7 — ! r E [ E (-l) lmAKl C R \N\(RUK))]x 

{ q ) KcN RcN\K 
(p) i (n-p-q) 

[ e (-lj^^^^U^h^ E l**X**l (DA) 

SCN\K \ q ) KCN 

The operator on the right-hand side of this equation is obviously positive 
definite. To complete the proof, we should check the normalization property 
that can be easily established by direct calculation ■ 
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Let us recast the last equation in a more convenient form that makes its 
structure more transparent. We have 

| i (p+q) 

A(n,p,q)^cO-«\*){*\ = 7^vE d q (Z) iP.h) 

P- { g ) ZCN 



where 
and 



d q (Z) = \* z )(* z \ (D.Q) 
(q) 

\*z) = E {-lf ZXm ^C z \ R \R) (D.7) 

RcZ 

Eqs.(D.5)-(D.7) may serve as a base for further analysis of both pure and en- 
semble representability problems. For example, at first glance the necessary 
condition of the ensemble representability of t q 

(ER) There exists an expansion ofd q = A(n,p, q)t q of the type ofEq.(D.5) 
such that for any Z C N,\Z\ = p + q the operator d q (Z) is nonnegative 

looks much stronger than simple positive semidefinitness of d q . 
The necessary and sufficient conditions of pure representability may be 
formulated with the aid of Eqs.(D.5)-(D.7). Let us introduce the set 

B n , p , q = {(Z, R) C N x N : \Z\ = p + qk\R\ = qkR C Z} (D.8) 

and the equivalence relation on this set 

(Z, R) ~ (Z', R') Z\R = Z'\R' (D.9) 

The set of all equivalence classes -B niPi9 contains (^j elements and in each 

equivalence class (Z, R) there are ( n ~ p ) elements. 

Theorem, q-electron operator t q is representable by pure p-electron state 
if and only if its image d q with respect to A(n,p,q) satisfies the following 
conditions: 

(i) There exists expansion 

i (p+q) 

dq = T^V E dq( Z ) (D.10) 



(V) 



ZCN 
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such that for each Z C N d q (Z) is either the null operator or corresponds to 
(unnormalized) pure q-electron state \^z); 
(ii) The mapping 

(Z,R) -> (-l)K z \ fl > nA "\*l(* z |i2) (D.ll) 

is constant on the equivalence classes (Z, R). 

If the conditions of this theorem are fullfilled then (up to normalization) 
the required pure state may be presented as 

|tf)= J2(-l)^ z \ R)nAN \^(^ z \R)\Z\R) (D.12) 

where the sum runs over equivalence classes. Unfortunately, condition (i) 
of the above theorem is depressingly non-constructive. Only for p + q = n 
conditions (i) and (ii) are trivial and the following corollary of the main 
theorem may be formulated. 

Corollary. For p + q = n q-electron operator t q is representable by 
pure p-electron state if and only if its image d q with respect to A(n,p,q) 
corresponds to pure q-electron state. 
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Appendix E. 

Let us show how by careful choice of the phase prefactors in Eq.(43) it is 
possible to ensure the unitary invariance of the convex set V n , P ,i- We have 



A(n,p,l)ef = (-!)""— "1*4 



n — p 



P jj 



n — p 

For arbitrary unitary operator u G U{J- n j) 

A(n, p, l)u ® «t e « = (-l)-*i-L_ ^ e «(i _ 

n P An 

^- e (-ir^ej^ fclfe < fe 



(£.1) 

(£.2) 



u ($u 



n ~P 

U(n,p,l)ef = (-ir -^-E4!(l-P^<J- 

' 4 P fci,fc 2 

(*l^fc2) 

A(n,p,l)«®«t e « = _(_i)«H^_^ e M 1ifci . u ^. + 

n P fci 



(£.3) 



(£.4) 



P 



E (-i) Qfeifc2 4 lfc2 ^^L- 



(£.5) 



u <g> v?A(n,p, l)e0 = (— 1) 
+(-l) ai 



: P 00 * , 

11 P fei 



2^ e u kli u k2j 



n-p - M 

(fei^fc 2 ) 



(£.6) 



From these equalities it is readily follows that to ensure the unitary invariance 
of A(n,p, 1) it is sufficient to require that either 



Q!00 = 0(mod 2) =>- otij = l(mod 2) 



(E.7a) 
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or 

o;00 = \{mod 2) =>- aij = 0(mod 2) (E.7b) 

for every i,j — l,...,n(i^ j). Our choice 0:00 = and = A^-jl 
for every 2 7^ j corresponds to the implication (E.7a). 
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Appendix F. 

The commutation relation (54) in coordinate form looks like 

(< E ] A% K {n,p,q-l) = A^ L (n,p,q) (F.l) 

KcL L'DK' 

where \L\ = \L'\ = q - s, and L, V C N\(I U J). 
Direct calculations show that 

(q-s-l) (p+\KnK'\-q\ (p+t-q\ , w v 

\- / u \KnK'\ \ \KnK'\ ) ( u t { t ){p-q){q-s) 

KCL \\KnK>\) { t ) P + t q 

and 

(q-s) (p+\LnL'\-q-l\ (p+t-q\ , ■ 

V ( U\LnL'\ { \Lnv\ ) _ / , t J (p-g)(n-p-g+l) 

h> ] " ( 1 (v) " 



(F.3) 



where t = \K' fl L|. 

Taking into account binomial prefactors in expression (43) leads readily 
to Eq.(F.l) that is equivalent to the commutation relation (54). 

It may be expected that there exist many reccurence relations of the type 
of Eq.(F.l) involving matrix elements of the operator A(n,p,q). 
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Appendix G. 

Lemma. Vectors w™ p (K) defined by Eq.(33) are lineary independent. 
Proof. Let us consider the vector equation 

E »kHU k ) = 

KCN\(IUJ) 

Using definitions (32) and (33), we can rewrite this equation in the form 

E ei J Tr[(» IJ )W p J lq (L)) = 

LcN\(IUJ) 

that is equivalent to the system of (™~ 2 /) of scalar equations 

Tri(f, Ij yw I p J lq (L)} = 0,LcN\(IUJ) 

But the set of vectors w p f q (L) is complete in the subspace of the operator 
space T n ^ q ® J 7 * determined by disjoint index sets / and J (see fll4"l)B 
Following general scheme of Sec. 4, we calculate 



<!■ „ P -g \ n P)- u p-g c u _ \q-sj\ p-g J 



H - -u 



p\ (n-q)\ 



E e^[G /J (n,p,g)]^, 

K'CN\(IUJ) 

(G.l) 
where 

Vg-s/ v v-i ) 

and s = |/| = |J|. Note that concrete blocks G IJ (n,p,q) thus introduced 
differ from blocks that could be constructed on the base of general definition 
(58) by non-essential combinatorial prefactor. The full-size block-diagonal 
matrix G(n,p,q) is 

q (s) fp-s\fn-q-s\ 

G(n, M )=9® ^ter C fJ (fl, P, 9) (G.3) 



s=0 7, J 



/ n-q 
p-q 



Since for any fixed I, J all blocks G IJ (n,p,q) are of similar structure we 
can confine ourselves to analysis of the "diagonal" case / = J — that 
corresponds to s = in Eqs.(G.l)-(G.3). 
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The problem of G ®(n,p, q) inversion can be reformulated as a combina- 
torial problem if we suppose that [G m (n,p, q)\~^ K u = f(n,p, q; \K fl K"\). It 
turns out, however, that we do not need [G(n, p, g)] -1 in itself but only the 
sums of the type 

b(n,p,q;\K'nL\) = ]T [G IJ (n,p,q)} K ) K (G.4) 

KcL 

as can be easily demonstrated by calculation of the expansion operator action 
on [G(n,p, q)]~H q . After rather complicated combinatorial manipulations we 
arrive at the linear system 

,£(;: 2 2 :;::)(;^;::-",)(n)(rJ(*'-'-^ 

(::;>- ^ 

where u = \K' fl L\. 

Its solution for all n,p,q (q < p < n) is 

b(n,p,q;r 2 ) = (-ir^h^^. (G.6) 

^ ' \q-ri) 

and our final formula necessary for explicit construction of section (60) is 

( v (p-s-\K'nL\-l\ 

e [G ij (n, P ,q)}-A = (-i) wi " ,nL| r_; V 9 T'r: L \ • ( G - 7 ) 
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Table 1: Results of simple calculations testing properties of the convex set 
V 2m , p ,2 : atoms 







Distance function 




Atomic 


Total 


deviation from the unit 


Absolute 


species 


energy 


The set 


The set 


error 




(a. u.) 




A- 1 (2m,p, 2)S 2m , 2 


in energy 
(a. u.) 


Li( 2 S) 


-7.433465 


< 10" 15 


< 2 • 10~ 8 


< 10~ 15 


Be( 1 5) 
B( 2 A) 


-14.618569 


4- 10" 13 


7- 10~ 8 


4-10" 12 


-24.591900 


4-10" 13 


4- 10- 8 


7- 10- 12 


C( 3 P) 


-37.761693 


< 10" 15 


< io~ 15 


< 10~ 15 


C^D) 
N( 4 S) 
0( 3 P) 


-37.707142 


< 10" 15 


6 • 10~ 15 


< io~ 15 


-54.480030 


< 10" 15 


1 • 10~ 15 


< io~ 15 


-74.844711 


< 10" 15 


< 10~ 15 


< 10~ 15 


OCD) 


-74.765527 


< 10" 15 


3 • 10" 15 


< io~ 15 
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Table 2: Results of simple calculations testing properties of the convex set 
V 2mjPi2 : small molecules and ions 



Molecular 
species 



Total 
energy 
(a. u.) 



Distance function 
deviation from the unit 



The set 



6- 10" 15 

4- lir 8 

2- 10- 7 
6- 10- 9 

4- lir 8 



The set 

A- 1 (2m,p,2)£ 2m:2 



KT 6 

io- 15 
io- 15 
< lir 15 
io- 14 
io- 15 
io- 14 
io- 14 
io- 15 
io- 15 



Absolute 

error 
in energy 
(a. u.) 



4-10- 14 
6 • 10~ 14 

1 ■ io- 13 

2 • 10~ 14 

3 • KT 13 

3 • 10- 13 

4 • 10- 13 
4 • KT 13 

1 • io- 13 

1 • KT 13 



LiH( 1 S+) 
CH 2 ( 3 5 2 ) 
CH 2 ( 1 A 1 ) 
NH 2 ( 2 5 2 ) 
NH 3 ( 1 A 1 ) 
H 2 0( 1 A 1 ) 
W( 3 B 2 ) 

H 2 0+( 2 A 1 ) 
H 2 0+( 2 5 2 ) 



-8.016132 

-38.979862 

-38.922201 

-55.478526 

-56.250080 

-76.077354 

-75.796065 

-75.765581 

-75.634415 

-75.560706 



1 



10" 
10" 



3- IO" 7 
6 • 10~ 10 
5- 10" 11 



5 
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